perm filename SUBSLM.F4[IRC,LCS] blob
sn#496778 filedate 1980-02-02 generic text, type T, neo UTF8
SUBROUTINE SS
COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
COMMON/SLM/A(2,200),B(2,200),C(2,200),I(200),I1(200),
1 I2(200),KT,TT
DATA L1/1/,L2/2/,L0/0/
C FIND MAX AND MIN POINTS AND SET SLOPES
ICOUNT=0
710 DO 100 K=2,N-1
KP=K-1
KQ=K+1
IF(Y(K).GT.Y(KP).AND.Y(K).GT.Y(KQ)) GO TO 20
IF(Y(K).LT.Y(KP).AND.Y(K).LT.Y(KQ)) GO TO 20
IF(X(K).GT.X(KP).AND.X(K).GT.X(KQ)) GO TO 30
IF(X(K).LT.X(KP).AND.X(K).LT.X(KQ)) GO TO 30
GO TO 100
20 S(K)=0
ICOUNT=ICOUNT+1
I(ICOUNT)=K
GO TO 100
30 S(K)=1001
ICOUNT=ICOUNT+1
I(ICOUNT)=K
100 CONTINUE
ICOUNT=ICOUNT+1
I(ICOUNT)=N
C FIND POINTS SUCH THAT X(K)=X(K+1) OR Y(K)=Y(K+1)
DO 900 K=2,N-1
IF(Y(K).EQ.Y(K+1)) I2(K)=1
IF(X(K).EQ.X(K+1)) I2(K)=1
900 CONTINUE
C SET THE SLOPES BETWEEN
IA=1
IX=0
110 IF(IX.EQ.ICOUNT) GO TO 200
IX=IX+1
IA=IA+1
120 IF(IA.EQ.I(IX)) GO TO 110
KP=IA+1
S(IA)=(Y(KP)-Y(IA-1))/(X(KP)-X(IA-1))
IA=KP
GO TO 120
200 CONTINUE
C SET BEGIN AND END SLOPES
CC SZ=1002001.
K=1
CALL S101
CALL S102(X(1),Y(1),X(2),Y(2),S(1),S(2))
IF(S(K).EQ.-1001.) S(K)=1001.
IF(S(K+1).EQ.-1001.) S(K+1)=1001.
K=N-1
CALL S101
CALL S102(X(N),Y(N),X(N-1),Y(N-1),S(N),S(N-1))
IF(S(K).EQ.-1001.) S(K)=1001.
IF(S(K+1).EQ.-1001.) S(K+1)=1001.
C RESET THE SLOPES FOR STRAIGHT LINES
DO 610 K=1,N-2
U1=SLM2(L0)
U2=SLM2(L1)
IF(ABS(U1-U2).GT..0001) GO TO 610
I1(K)=1
I1(K+1)=1
S(K)=U1
S(K+1)=U1
S(K+2)=U1
610 CONTINUE
C ADD POINTS
K1=N-1
DO 300 K=1,K1
IF(I1(K).EQ.1) GO TO 300
IF(I2(K).EQ.1) GO TO 300
FLG=0.
IF(S(K).EQ.0..AND.S(K+1).EQ.0.) FLG=1.
IF(S(K).EQ.1001..AND.S(K+1).EQ.1001.) FLG=2.
IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 205
CC SZ=1002001.
CALL S101
KP=K+1
U=(Y(K)-Y(KP)-S(K)*(X(K)-X(KP)))/(S(KP)-S(K))
AX=X(KP)+U
AY=Y(KP)+S(KP)*U
XA=X(K)
1610 XB=X(KP)
IF(XA.LE.XB) GO TO 202
XA=X(KP)
XB=X(K)
202 YA=Y(K)
YB=Y(KP)
IF(YA.LE.YB) GO TO 204
YA=Y(KP)
YB=Y(K)
204 CONTINUE
IF(AX.GE.XA.AND.AX.LE.XB.AND.
1 AY.GE.YA.AND.AY.LE.YB) GO TO 290
205 K1=K1+1
DO 210 K2=K+1,N
K3=N+K+1-K2
KP=K3+1
X(KP)=X(K3)
Y(KP)=Y(K3)
I1(KP)=I1(K3)
I2(KP)=I2(K3)
210 S(KP)=S(K3)
IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 280
N=N+1
KP=K+1
KQ=K+2
XK=X(K)
XQ=X(KQ)
XZ=XK-XQ
C Z0=X(K)**2-X(KQ)**2
Z0=XK*XK+XQ*XQ
C Z1=Y(K)-Y(KQ)-S(K)*(X(K)-X(KQ))-((X(K)+X(KQ))/2.
C 1 -X(K))*(S(K)-S(KQ))
Z1=Y(K)-Y(KQ)-S(K)*XZ-(XK+XQ)/2.-XK*(S(K)-S(KQ))
C Z2=X(K)**3-X(KQ)**3-3*(X(K)-X(KQ))*X(K)**2
C 1 -1.5*(X(K)+X(KQ))*Z0+3*X(K)*Z0
Z2=XK**3-XQ**3-3*(XK-XQ)*XK*XK-1.5*(XK+XQ)*Z0+3*XK*Z0
AZ=Z1/Z2
C BZ=(S(K)-S(KQ)-3*(X(K)**2-X(KQ)**2)*AZ)/
C 1 (2*(X(K)-X(KQ)))
BZ=(S(K)-S(KQ)-3*(XK*XK-XQ*XQ)*AZ)/(2*XZ)
C CZ=S(K)-3*AZ*X(K)**2-2*BZ*X(K)
CZ=S(K)-3*AZ*XK*XK-2*BZ*XK
C DZ=Y(K)-AZ*X(K)**3-BZ*X(K)**2-CZ*X(K)
DZ=Y(K)-AZ*XK**3-BZ*XK*XK-CZ*XK
X(KP)=-BZ/(3.*AZ)
Y(KP)=AZ*X(KP)**3+BZ*X(KP)**2+CZ*X(KP)+DZ
S(KP)=3*AZ*X(KP)**2+2*BZ*X(KP)+CZ
IF(Y(KP).LT.YB.AND.Y(KP).GT.YA) GO TO 290
X(KP)=(XQ+XK)/2.
Y(KP)=(Y(KQ)+Y(K))/2.
S(KP)=0.
GO TO 290
280 N=N+1
IF(FLG.EQ.2.) GO TO 282
S(K+1)=2.*(Y(K+2)-Y(K))/(X(K+2)-X(K))
GO TO 284
282 S(K+1)=(Y(K+2)-Y(K))/(2.*(X(K+2)-X(K)))
284 X(K+1)=(X(K+2)+X(K))/2.
Y(K+1)=(Y(K+2)+Y(K))/2.
290 IF(S(K).EQ.-1001.) S(K)=1001.
IF(S(K+1).EQ.-1001.) S(K+1)=1001.
300 CONTINUE
C FIT THE POINTS WITH N-1 CURVES
305 DO 400 K=1,N-1
KP=K+1
IF(I1(K).EQ.0) GO TO 310
A(1,K)=0.
A(2,K)=0.
B(1,K)=X(KP)-X(K)
B(2,K)=Y(KP)-Y(K)
C(1,K)=X(K)
C(2,K)=Y(K)
GO TO 400
310 CALL S101
B(1,K)=2*(Y(KP)-Y(K)-S(KP)*(X(KP)-X(K)))/
1 (S(K)-S(KP))
B(2,K)=S(K)*B(1,K)
A(1,K)=X(KP)-X(K)-B(1,K)
A(2,K)=Y(KP)-Y(K)-S(K)*B(1,K)
C(1,K)=X(K)
C(2,K)=Y(K)
400 CONTINUE
C CALCULATE 512 POINTS
M=N-1
R=M/511.
T=-R
DO 500 L=1,511
T=T+R
KT=T+1
KU=T
TT=T-KU
X1(L)=SLM3(L1)
500 Y1(L)=SLM3(L2)
X1(512)=X(N)
Y1(512)=Y(N)
600 RETURN
END
C GIVEN TWO POINTS (XX1,YY1), (XX2,YY2) AND A SLOPE SS2
C THIS ROUTINE FINDS THE SLOPE SS1
SUBROUTINE S102(XX1,YY1,XX2,YY2,SS1,SS2)
LF=0
IF(SS2.NE.0.) GO TO 5
IF((XX1.LT.XX2.AND.YY1.LT.YY2).OR.
1 (XX2.LT.XX1.AND.YY2.LT.YY1)) SS1=1001.
IF((XX1.LT.XX2.AND.YY1.GT.YY2).OR.
1 (XX2.LT.XX1.AND.YY2.GT.YY1)) SS1=-1001.
GO TO 100
5 T=(YY1-YY2)/SS2
X=XX2+T
Y=YY1
XMAX=XX2
XMIN=XX1
IF(XX2.GE.XX1) GO TO 10
XMAX=XX1
XMIN=XX2
10 IF(X.GT.XMIN.AND.X.LT.XMAX) GO TO 20
T=XX1-XX2
Y=YY2+SS2*T
X=XX1
LF=1
20 IF(LF.EQ.1) GO TO 30
D2=XX1-X
D1=X-XX2
R=D2/D1
XX=XX2
YY=(YY2+R*YY1)/(1.+R)
GO TO 40
30 D2=YY1-Y
D1=Y-YY2
R=D2/D1
YY=YY2
XX=(XX2+R*XX1)/(1.+R)
40 SS1=(YY-YY1)/(XX-XX1)
100 RETURN
END
SUBROUTINE S101
COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
DATA SZ/1002001./
SS=S(K)**2
SSS=S(K+1)**2
IF(SS.NE.SZ.AND.SSS.NE.SZ)RETURN
IF(SS.EQ.SZ.AND.X(K).LT.X(K+1).AND.
1 Y(K).GT.Y(K+1)) S(K)=-1001.
IF(S(K)**2.EQ.SZ.AND.X(K).GT.X(K+1).AND.
1 Y(K+1).GT.Y(K)) S(K)=-1001.
IF(SSS.EQ.SZ.AND.X(K).LT.X(K+1).AND.
1 Y(K).GT.Y(K+1)) S(K+1)=-1001.
IF(S(K+1)**2.EQ.SZ.AND.X(K).GT.X(K+1).AND.
1 Y(K).LT.Y(K+1)) S(K+1)=-1001.
100 RETURN
END
FUNCTION SLM2(L)
COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
M=K+L
J=M+1
SLM2=(Y(J)-Y(M))/(X(J)-X(M))
END
FUNCTION SLM3(L)
COMMON/SLM/A(2,200),B(2,200),C(2,200),I(200),I1(200),
1 I2(200),KT,TT
SLM3=A(L,KT)*TT**2+B(L,KT)*TT+C(L,KT)
END